Analysis of Freesurfer parcellation Genetic Template, surface area and cortical thickness differences

library(easypackages)
libraries("here","heplots","ggplot2","ggseg","ggsegExtra","ggsegChen","tidyverse","patchwork","pheatmap","psych")
source(here("code","cohens_d.R"))
source(here("code","get_ggColorHue.R"))

datapath = here("data","tidy")
plotdir = here("plots")

pheno_data_gex =read.csv(file.path(datapath,"labelData_allGEXMRIsubs.csv"))
pheno_data_mri =read.csv(file.path(datapath,"labelData_all_MRI.csv"))

# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)

fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)

# region_names = colnames(fs_data)[5:ncol(fs_data)]
# region_names = c("CortexVol","SAtotal","CTmean",region_names)
region_names = c("CortexVol","SAtotal","CTmean")

fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")


#------------------------------------------------------------------------------
# Remove the longitudinal duplicates
dup_df = data.frame(t(table(fs_data$subjectId)))
colnames(dup_df)[2:3] = c("subjectId","n")
dup_subs = dup_df$subjectId[dup_df$n>1]

rows2keep = numeric()
rows2remove = numeric()
for (isub in dup_subs){
  dup_indx = which(fs_data$subjectId==isub)
  rows2keep = c(rows2keep,dup_indx[1])
  rows2remove = c(rows2remove,dup_indx[2:length(dup_indx)])
}

fs_data = fs_data[-rows2remove,] 
#------------------------------------------------------------------------------


new_fs_data = fs_data

colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
                 "tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
                 "pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use

for (icol in 1:length(region_names)){
  
  y_var = region_names[icol]
  
  if(is.element(y_var,c("CortexVol","SAtotal","CTmean"))){
    # make formula
    form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
    full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
      covname2use = c("sexMale","scan_age")
  }else{
    # make formula
    form2use = as.formula(sprintf("%s ~ Group + CortexVol + sex + scan_age", y_var))
    full_model = model.matrix(~0+ Group + CortexVol + sex + scan_age, data = fs_data)
    covname2use = c("CortexVol","sexMale","scan_age")
    # #--------------------------------------------------------------------------
    # # remove CortexVol as covariate since GCLUST is already globally controlled for
    # # make formula
    # form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
    # full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
    # covname2use = c("sexMale","scan_age")
    # #--------------------------------------------------------------------------
  }
  
  # run linear model
  mod2use = lm(formula = form2use, data = fs_data)
  
  # run an ANOVA
  res2use = anova(mod2use)
  
  # grab ANOVA results and store into output_res
  output_res[y_var, "Fstat"] = res2use["Group", "F value"]
  output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
  
  # compute partial eta-squared as the effect size of the interaction effect
  eta_sq_res = etasq(mod2use)
  output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
  
  
    # remove covariates before effect size computation
    beta1 = mod2use$coefficients[covname2use, drop = FALSE]
    beta1[is.na(beta1)] = 0

    new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
    
    # compute effect sizes
    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"] 
    
    mask1 = fs_data$Group=="Poor"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 


    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="Poor"
    output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])

    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 

}

output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
output_res[order(output_res$pval),]
##               Fstat         pval       etasq         fdrq es.TD_vs_Good
## SAtotal   15.397759 1.143503e-06 0.072173002 3.430509e-06   -0.27132289
## CortexVol 14.303102 2.741889e-06 0.075345308 4.112833e-06   -0.26212136
## CTmean     2.803486 6.464027e-02 0.002252612 6.464027e-02    0.07724623
##           es.TD_vs_Poor es.Good_vs_Poor tstat.TD_vs_Good tstat.TD_vs_Poor
## SAtotal      -0.7113176     -0.40788261       -1.2806955       -2.8440667
## CortexVol    -0.7087278     -0.42651049       -1.1088191       -2.9849992
## CTmean        0.1155233      0.03662885        0.3724289        0.1765433
##           tstat.Good_vs_Poor pval.TD_vs_Good pval.TD_vs_Poor pval.Good_vs_Poor
## SAtotal            1.7936785       0.2039557     0.005637100        0.07706336
## CortexVol          1.8842593       0.2707896     0.003748508        0.06356704
## CTmean            -0.2434033       0.7105462     0.860307799        0.80838511
##           fdrq.TD_vs_Good fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## SAtotal         0.4061844     0.008455649         0.1155950
## CortexVol       0.4061844     0.008455649         0.1155950
## CTmean          0.7105462     0.860307799         0.8083851

Make plots of total cortical volume, surface area, and mean cortical thickness

dotSize = 6
fontSize=28

fdr_thresh = 0.05
output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
output_res_fdr[order(output_res_fdr$pval),]
##              Fstat         pval      etasq         fdrq es.TD_vs_Good
## SAtotal   15.39776 1.143503e-06 0.07217300 3.430509e-06    -0.2713229
## CortexVol 14.30310 2.741889e-06 0.07534531 4.112833e-06    -0.2621214
##           es.TD_vs_Poor es.Good_vs_Poor tstat.TD_vs_Good tstat.TD_vs_Poor
## SAtotal      -0.7113176      -0.4078826        -1.280695        -2.844067
## CortexVol    -0.7087278      -0.4265105        -1.108819        -2.984999
##           tstat.Good_vs_Poor pval.TD_vs_Good pval.TD_vs_Poor pval.Good_vs_Poor
## SAtotal             1.793679       0.2039557     0.005637100        0.07706336
## CortexVol           1.884259       0.2707896     0.003748508        0.06356704
##           fdrq.TD_vs_Good fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## SAtotal         0.4061844     0.008455649          0.115595
## CortexVol       0.4061844     0.008455649          0.115595
# Plot Total Cortical Volume
cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CortexVol"
reg_name = "Total Cortical Volume"
metric_name = "Volume"
ylims = c(450000,700000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = here("plots","TotalCorticalVolume.pdf"))
p1

# Plot Total Surface Area
cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "SAtotal"
reg_name = "Total Surface Area"
metric_name = "SA"
ylims = c(90000,220000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = here("plots","TotalSurfaceArea.pdf"))
p1

# Plot Mean Cortical Thickness
cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CTmean"
reg_name = "Mean Cortical Thickness"
metric_name = "CT"
ylims = c(2.75,3.75)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = here("plots","MeanCorticalThickness.pdf"),width=7,height=5)
p1

Do modeling of local differences without adjustment for total cortical volume because GCLUST is already globally corrected.

# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)

fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)

region_names = colnames(fs_data)[5:ncol(fs_data)]
# region_names = c("CortexVol","SAtotal","CTmean",region_names)

fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")

#------------------------------------------------------------------------------
# Remove duplicates
dup_df = data.frame(t(table(fs_data$subjectId)))
colnames(dup_df)[2:3] = c("subjectId","n")
dup_subs = dup_df$subjectId[dup_df$n>1]

rows2keep = numeric()
rows2remove = numeric()
for (isub in dup_subs){
  dup_indx = which(fs_data$subjectId==isub)
  rows2keep = c(rows2keep,dup_indx[1])
  rows2remove = c(rows2remove,dup_indx[2:length(dup_indx)])
}

fs_data = fs_data[-rows2remove,] 
#------------------------------------------------------------------------------


new_fs_data = fs_data

colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
                 "tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
                 "pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use

for (icol in 1:length(region_names)){
  
  y_var = region_names[icol]
  
  # make formula
  form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
  full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
  covname2use = c("sexMale","scan_age")

  # run linear model
  mod2use = lm(formula = form2use, data = fs_data)
  
  # run an ANOVA
  res2use = anova(mod2use)
  
  # grab ANOVA results and store into output_res
  output_res[y_var, "Fstat"] = res2use["Group", "F value"]
  output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
  
  # compute partial eta-squared as the effect size of the interaction effect
  eta_sq_res = etasq(mod2use)
  output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
  
  
    # remove covariates before effect size computation
    beta1 = mod2use$coefficients[covname2use, drop = FALSE]
    beta1[is.na(beta1)] = 0

    new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
    
    # compute effect sizes
    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"] 
    
    mask1 = fs_data$Group=="Poor"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 


    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="Poor"
    output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])

    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 

}

output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
write.csv(output_res, file = here("results","GCLUST_regional_group_differences_stats.csv"))
output_res[order(output_res$pval),]
##                                         Fstat        pval        etasq
## middletemporal.rh.thickness        5.94103928 0.003480438 0.0734094108
## posterolateraltemporal.rh.area     5.46548677 0.005369823 0.0598876648
## anteromedialtemporal.lh.area       5.15535591 0.007137177 0.0658046617
## superiorparietal.rh.area           4.07830238 0.019379782 0.0287779414
## medialprefrontal.lh.thickness      4.06769013 0.019573122 0.0275112055
## middletemporal.lh.thickness        3.41762916 0.036069477 0.0343140556
## superiorparietal.lh.area           3.31498536 0.039747640 0.0330204534
## dorsomedialfrontal.lh.area         3.15695112 0.046171733 0.0404887240
## medialprefrontal.rh.thickness      2.92870020 0.057364440 0.0129543380
## temporalpole.rh.thickness          2.76347611 0.067157819 0.0229000518
## orbitofrontal.rh.area              2.56863271 0.080920331 0.0853698327
## temporalpole.lh.thickness          2.38151980 0.096838507 0.0370636668
## occipital.lh.area                  2.24561269 0.110368376 0.0271548766
## posterolateraltemporal.lh.area     1.93258720 0.149327978 0.0140083037
## inferiorparietal.lh.area           1.89942858 0.154201588 0.0328601053
## dorsolateralprefrotal.lh.thickness 1.86898803 0.158818058 0.0112166280
## inferiorparietal.lh.thickness      1.84940380 0.161862175 0.0293629309
## orbitofrontal.lh.area              1.77668475 0.173694196 0.0542587241
## perisylvian.lh.thickness           1.76919001 0.174962619 0.0083808211
## superiortemporal.rh.area           1.70160553 0.186834976 0.0080576156
## dorsolateralprefrotal.rh.thickness 1.69970989 0.187179541 0.0113713568
## anteromedialtemporal.rh.area       1.68188886 0.190450595 0.0078339922
## dorsomedialfrontal.rh.area         1.47819586 0.232240462 0.0313502126
## inferiorparietal.rh.area           1.46476311 0.235304210 0.0322111775
## motor_premotor_SMA.rh.thickness    1.30638497 0.274684555 0.0447032529
## motor_premotor_SMA.lh.thickness    1.18615421 0.309007959 0.0610710194
## perisylvian.rh.thickness           1.15784633 0.317705086 0.0132777649
## ventralfrontal.rh.thickness        1.12958718 0.326635646 0.0037612076
## ventromedialoccipital.lh.thickness 1.03790036 0.357408138 0.0338876279
## inferiorparietal.rh.thickness      1.03713528 0.357676959 0.0196799443
## ventromedialoccipital.rh.thickness 0.87775545 0.418411359 0.0162415410
## superiorparietal.rh.thickness      0.63051356 0.534102064 0.0105049232
## medialtemporal.lh.thickness        0.59903257 0.551004680 0.0101886096
## dorsolateralprefrontal.lh.area     0.54820187 0.579453005 0.0117343658
## motorpremotor.lh.area              0.50722783 0.603469423 0.0069401882
## superiorparietal.lh.thickness      0.47653492 0.622120938 0.0080707604
## dorsolateralprefrontal.rh.area     0.43942349 0.645457966 0.0040848178
## occipital.rh.area                  0.37917411 0.685257301 0.0066204803
## occipital.rh.thickness             0.32871018 0.720508629 0.0041295970
## motorpremotor.rh.area              0.32742064 0.721433208 0.0105939533
## ventralfrontal.lh.thickness        0.30587253 0.737062641 0.0030880881
## precuneus.rh.area                  0.28560197 0.752079731 0.0035408362
## occipital.lh.thickness             0.28391299 0.753344940 0.0054019402
## medialtemporal.rh.thickness        0.21417637 0.807519105 0.0003133628
## parsopercularis.lh.area            0.15301539 0.858286507 0.0053473121
## parsopercularis.rh.area            0.15014302 0.860749032 0.0008351415
## precuneus.lh.area                  0.08999012 0.914002879 0.0017483701
## superiortemporal.lh.area           0.01223912 0.987836730 0.0029551953
##                                         fdrq es.TD_vs_Good es.TD_vs_Poor
## middletemporal.rh.thickness        0.1141948  -0.663418297  -0.392605943
## posterolateraltemporal.rh.area     0.1141948   0.030295479  -0.535892781
## anteromedialtemporal.lh.area       0.1141948   0.105637215  -0.533075090
## superiorparietal.rh.area           0.1879020  -0.101374144   0.311313007
## medialprefrontal.lh.thickness      0.1879020   0.118687302   0.420863328
## middletemporal.lh.thickness        0.2725552  -0.149600918  -0.519216132
## superiorparietal.lh.area           0.2725552  -0.141813175   0.330575141
## dorsomedialfrontal.lh.area         0.2770304  -0.360530708   0.121131105
## medialprefrontal.rh.thickness      0.3059437   0.108139104   0.282925991
## temporalpole.rh.thickness          0.3223575  -0.294484547  -0.368873339
## orbitofrontal.rh.area              0.3531069   0.669848517   0.631671827
## temporalpole.lh.thickness          0.3873540   0.022831450  -0.431998283
## occipital.lh.area                  0.4075140   0.417678877   0.183761925
## posterolateraltemporal.lh.area     0.4155286  -0.108519406  -0.306559251
## inferiorparietal.lh.area           0.4155286  -0.215460672  -0.459224327
## dorsolateralprefrotal.lh.thickness 0.4155286   0.119445339   0.275104346
## inferiorparietal.lh.thickness      0.4155286   0.188388475  -0.258606627
## orbitofrontal.lh.area              0.4155286   0.214206795   0.585759310
## perisylvian.lh.thickness           0.4155286  -0.210011882  -0.170481439
## superiortemporal.rh.area           0.4155286   0.209040493   0.153640728
## dorsolateralprefrotal.rh.thickness 0.4155286   0.215372099  -0.006725081
## anteromedialtemporal.rh.area       0.4155286  -0.096075705  -0.232557068
## dorsomedialfrontal.rh.area         0.4706084   0.060393743   0.482914856
## inferiorparietal.rh.area           0.4706084  -0.464685948  -0.104099702
## motor_premotor_SMA.rh.thickness    0.5273943   0.319915972   0.547023221
## motor_premotor_SMA.lh.thickness    0.5599468   0.448981145   0.648166073
## perisylvian.rh.thickness           0.5599468   0.078911716  -0.206641773
## ventralfrontal.rh.thickness        0.5599468  -0.148932409  -0.063507467
## ventromedialoccipital.lh.thickness 0.5722831  -0.304107726  -0.503207370
## inferiorparietal.rh.thickness      0.5722831   0.330777772   0.055870712
## ventromedialoccipital.rh.thickness 0.6478627  -0.293894577  -0.261454543
## superiorparietal.rh.thickness      0.8011531   0.219385212   0.012428666
## medialtemporal.lh.thickness        0.8014614   0.044354354   0.241938279
## dorsolateralprefrontal.lh.area     0.8180513   0.122140978   0.270635677
## motorpremotor.lh.area              0.8276152  -0.032713682   0.165653399
## superiorparietal.lh.thickness      0.8294946  -0.100464989   0.138091765
## dorsolateralprefrontal.rh.area     0.8373509   0.005594604  -0.132562680
## occipital.rh.area                  0.8409432   0.157964681   0.183810743
## occipital.rh.thickness             0.8409432  -0.012097343   0.140113306
## motorpremotor.rh.area              0.8409432  -0.232932678  -0.188478535
## ventralfrontal.lh.thickness        0.8409432  -0.135651383  -0.111708180
## precuneus.rh.area                  0.8409432  -0.015804628   0.118539861
## occipital.lh.thickness             0.8409432  -0.061960608   0.118584159
## medialtemporal.rh.thickness        0.8809299  -0.033634161   0.006069805
## parsopercularis.lh.area            0.8981729  -0.049461374  -0.206013379
## parsopercularis.rh.area            0.8981729  -0.057073500   0.001465130
## precuneus.lh.area                  0.9334497  -0.026753735   0.075252332
## superiortemporal.lh.area           0.9878367  -0.087919189  -0.137690521
##                                    es.Good_vs_Poor tstat.TD_vs_Good
## middletemporal.rh.thickness            0.308293704      -2.85110017
## posterolateraltemporal.rh.area        -0.578778313       0.36311947
## anteromedialtemporal.lh.area          -0.627195070       0.37894890
## superiorparietal.rh.area               0.438346980      -0.63925964
## medialprefrontal.lh.thickness          0.274891792       0.41679478
## middletemporal.lh.thickness           -0.296831346      -0.66119236
## superiorparietal.lh.area               0.440699732      -0.67574253
## dorsomedialfrontal.lh.area             0.480764273      -1.63798702
## medialprefrontal.rh.thickness          0.175740176       0.64582311
## temporalpole.rh.thickness             -0.057108681      -1.34457649
## orbitofrontal.rh.area                 -0.047911120       2.95765605
## temporalpole.lh.thickness             -0.430949000       0.05637093
## occipital.lh.area                     -0.221406959       1.73969339
## posterolateraltemporal.lh.area        -0.187309136      -0.49040595
## inferiorparietal.lh.area              -0.242790899      -0.91920137
## dorsolateralprefrotal.lh.thickness     0.144043968       0.39221574
## inferiorparietal.lh.thickness         -0.426697227       1.03758906
## orbitofrontal.lh.area                  0.385080584       0.92711527
## perisylvian.lh.thickness               0.035500517      -0.94804045
## superiortemporal.rh.area              -0.056252035       1.09667446
## dorsolateralprefrotal.rh.thickness    -0.245041864       0.62950223
## anteromedialtemporal.rh.area          -0.121524702      -0.54955871
## dorsomedialfrontal.rh.area             0.330214373       0.25946827
## inferiorparietal.rh.area               0.298412101      -2.24081659
## motor_premotor_SMA.rh.thickness        0.215161604       1.52075062
## motor_premotor_SMA.lh.thickness        0.169031261       2.12900175
## perisylvian.rh.thickness              -0.291323899       0.25516795
## ventralfrontal.rh.thickness            0.086835099      -0.89476216
## ventromedialoccipital.lh.thickness    -0.122955125      -1.32562092
## inferiorparietal.rh.thickness         -0.279450323       1.52372353
## ventromedialoccipital.rh.thickness     0.046647076      -1.22983796
## superiorparietal.rh.thickness         -0.230412846       1.11877053
## medialtemporal.lh.thickness            0.201395580       0.22010597
## dorsolateralprefrontal.lh.area         0.151199701       0.67279733
## motorpremotor.lh.area                  0.193878295      -0.20536507
## superiorparietal.lh.thickness          0.209895159      -0.28637751
## dorsolateralprefrontal.rh.area        -0.150628075       0.09706286
## occipital.rh.area                      0.046920394       0.56402834
## occipital.rh.thickness                 0.144938773      -0.04980532
## motorpremotor.rh.area                  0.049542202      -0.94441590
## ventralfrontal.lh.thickness            0.008053005      -0.82119741
## precuneus.rh.area                      0.158720852      -0.11229016
## occipital.lh.thickness                 0.191400339      -0.37207575
## medialtemporal.rh.thickness            0.039101046      -0.03002987
## parsopercularis.lh.area               -0.123165860      -0.29915811
## parsopercularis.rh.area                0.066652372      -0.21632630
## precuneus.lh.area                      0.104640544      -0.07671375
## superiortemporal.lh.area              -0.046795143      -0.06771132
##                                    tstat.TD_vs_Poor tstat.Good_vs_Poor
## middletemporal.rh.thickness             -1.55784151        -1.27463943
## posterolateraltemporal.rh.area          -2.18864474         2.70658431
## anteromedialtemporal.lh.area            -1.89885775         2.70178261
## superiorparietal.rh.area                 1.52415173        -1.92100876
## medialprefrontal.lh.thickness            1.79227030        -1.20242122
## middletemporal.lh.thickness             -1.94016360         1.36426487
## superiorparietal.lh.area                 1.51281550        -1.83973682
## dorsomedialfrontal.lh.area               0.41544508        -2.20890355
## medialprefrontal.rh.thickness            1.10578637        -0.69398289
## temporalpole.rh.thickness               -1.35147102         0.46949803
## orbitofrontal.rh.area                    2.90748713         0.25800500
## temporalpole.lh.thickness               -1.76002661         1.96294321
## occipital.lh.area                        0.52774695         0.85917149
## posterolateraltemporal.lh.area          -0.95004682         0.94694197
## inferiorparietal.lh.area                -1.96219360         1.06956480
## dorsolateralprefrotal.lh.thickness       1.35489026        -0.59374712
## inferiorparietal.lh.thickness           -1.18684886         1.84689550
## orbitofrontal.lh.area                    2.43578917        -1.69232724
## perisylvian.lh.thickness                -0.86889078        -0.18510829
## superiortemporal.rh.area                 0.32997478         0.25054848
## dorsolateralprefrotal.rh.thickness       0.14948667         0.95760379
## anteromedialtemporal.rh.area            -0.86457779         0.49201307
## dorsomedialfrontal.rh.area               2.12769481        -1.38047753
## inferiorparietal.rh.area                -0.16728797        -1.24213808
## motor_premotor_SMA.rh.thickness          2.23692644        -0.89401899
## motor_premotor_SMA.lh.thickness          2.80489444        -0.65068283
## perisylvian.rh.thickness                -0.94393787         1.15847060
## ventralfrontal.rh.thickness             -0.27982984        -0.45174958
## ventromedialoccipital.lh.thickness      -2.39680348         0.37093793
## inferiorparietal.rh.thickness           -0.20881417         0.99827580
## ventromedialoccipital.rh.thickness      -0.86594446        -0.15440560
## superiorparietal.rh.thickness            0.06444381         1.01733021
## medialtemporal.lh.thickness              1.07207894        -0.72919504
## dorsolateralprefrontal.lh.area           0.90469439        -0.65364085
## motorpremotor.lh.area                    0.88366726        -0.79375138
## superiorparietal.lh.thickness            0.57336220        -0.93346662
## dorsolateralprefrontal.rh.area          -0.45090801         0.68606122
## occipital.rh.area                        0.57617753        -0.31328688
## occipital.rh.thickness                   0.57615386        -0.67486217
## motorpremotor.rh.area                   -0.98407355        -0.24873178
## ventralfrontal.lh.thickness             -0.35130305        -0.03719715
## precuneus.rh.area                        0.54329017        -0.82107902
## occipital.lh.thickness                   0.25107936        -1.01417277
## medialtemporal.rh.thickness              0.24565262         0.06611091
## parsopercularis.lh.area                 -0.92713325         0.49447561
## parsopercularis.rh.area                 -0.31262077        -0.39751903
## precuneus.lh.area                        0.56040020        -0.42510178
## superiortemporal.lh.area                -1.14151848         0.17263804
##                                    pval.TD_vs_Good pval.TD_vs_Poor
## middletemporal.rh.thickness            0.005525174     0.123169627
## posterolateraltemporal.rh.area         0.717461549     0.031502420
## anteromedialtemporal.lh.area           0.705717336     0.061143805
## superiorparietal.rh.area               0.524457792     0.131365824
## medialprefrontal.lh.thickness          0.677931016     0.076823579
## middletemporal.lh.thickness            0.510365327     0.055839716
## superiorparietal.lh.area               0.501128867     0.134218641
## dorsomedialfrontal.lh.area             0.105303557     0.678914597
## medialprefrontal.rh.thickness          0.520219365     0.272092547
## temporalpole.rh.thickness              0.182515299     0.180307013
## orbitofrontal.rh.area                  0.004061395     0.004699057
## temporalpole.lh.thickness              0.955185137     0.082178529
## occipital.lh.area                      0.085711306     0.599117724
## posterolateraltemporal.lh.area         0.625172512     0.344914593
## inferiorparietal.lh.area               0.360719517     0.053174065
## dorsolateralprefrotal.lh.thickness     0.695928937     0.179219372
## inferiorparietal.lh.thickness          0.302548979     0.238756489
## orbitofrontal.lh.area                  0.356621001     0.017053979
## perisylvian.lh.thickness               0.345928755     0.387474182
## superiortemporal.rh.area               0.276033472     0.742270900
## dorsolateralprefrotal.rh.thickness     0.530791893     0.881541209
## anteromedialtemporal.rh.area           0.584134133     0.389823804
## dorsomedialfrontal.rh.area             0.795932285     0.036403210
## inferiorparietal.rh.area               0.027775608     0.867560414
## motor_premotor_SMA.rh.thickness        0.132216658     0.028039496
## motor_premotor_SMA.lh.thickness        0.036291546     0.006299476
## perisylvian.rh.thickness               0.799240172     0.348008513
## ventralfrontal.rh.thickness            0.373565377     0.780321412
## ventromedialoccipital.lh.thickness     0.188691907     0.018841738
## inferiorparietal.rh.thickness          0.131472705     0.835117540
## ventromedialoccipital.rh.thickness     0.222317948     0.389078318
## superiorparietal.rh.thickness          0.266544921     0.948775623
## medialtemporal.lh.thickness            0.826342416     0.286868405
## dorsolateralprefrontal.lh.area         0.502991170     0.368310330
## motorpremotor.lh.area                  0.837802144     0.379491172
## superiorparietal.lh.thickness          0.775320195     0.567987958
## dorsolateralprefrontal.rh.area         0.922916239     0.653260634
## occipital.rh.area                      0.574293259     0.566092762
## occipital.rh.thickness                 0.960400139     0.566108683
## motorpremotor.rh.area                  0.347765765     0.328009702
## ventralfrontal.lh.thickness            0.413944279     0.726273163
## precuneus.rh.area                      0.910871155     0.588422089
## occipital.lh.thickness                 0.710808117     0.802388604
## medialtemporal.rh.thickness            0.976117097     0.806572522
## parsopercularis.lh.area                0.765585637     0.356611722
## parsopercularis.rh.area                0.829277289     0.755372241
## precuneus.lh.area                      0.939040433     0.576753254
## superiortemporal.lh.area               0.946182305     0.257017722
##                                    pval.Good_vs_Poor fdrq.TD_vs_Good
## middletemporal.rh.thickness              0.206535221       0.1326042
## posterolateraltemporal.rh.area           0.008483109       0.9761171
## anteromedialtemporal.lh.area             0.008594918       0.9761171
## superiorparietal.rh.area                 0.058690208       0.9761171
## medialprefrontal.lh.thickness            0.233139903       0.9761171
## middletemporal.lh.thickness              0.176733527       0.9761171
## superiorparietal.lh.area                 0.069929008       0.9761171
## dorsomedialfrontal.lh.area               0.030364342       0.7933000
## medialprefrontal.rh.thickness            0.489926021       0.9761171
## temporalpole.rh.thickness                0.640132452       0.9057212
## orbitofrontal.rh.area                    0.797139090       0.1326042
## temporalpole.lh.thickness                0.053516450       0.9761171
## occipital.lh.area                        0.393097427       0.7933000
## posterolateraltemporal.lh.area           0.346835702       0.9761171
## inferiorparietal.lh.area                 0.288387881       0.9437441
## dorsolateralprefrotal.lh.thickness       0.554542121       0.9761171
## inferiorparietal.lh.thickness            0.068871440       0.9437441
## orbitofrontal.lh.area                    0.094907165       0.9437441
## perisylvian.lh.thickness                 0.853664547       0.9437441
## superiortemporal.rh.area                 0.802876807       0.9437441
## dorsolateralprefrotal.rh.thickness       0.341467101       0.9761171
## anteromedialtemporal.rh.area             0.624206920       0.9761171
## dorsomedialfrontal.rh.area               0.171710503       0.9761171
## inferiorparietal.rh.area                 0.218216386       0.4354985
## motor_premotor_SMA.rh.thickness          0.374290497       0.7933000
## motor_premotor_SMA.lh.thickness          0.517322378       0.4354985
## perisylvian.rh.thickness                 0.250500151       0.9761171
## ventralfrontal.rh.thickness              0.652807328       0.9437441
## ventromedialoccipital.lh.thickness       0.711772955       0.9057212
## inferiorparietal.rh.thickness            0.321488880       0.7933000
## ventromedialoccipital.rh.thickness       0.877721941       0.9437441
## superiorparietal.rh.thickness            0.312402770       0.9437441
## medialtemporal.lh.thickness              0.468247853       0.9761171
## dorsolateralprefrontal.lh.area           0.515425366       0.9761171
## motorpremotor.lh.area                    0.429948432       0.9761171
## superiorparietal.lh.thickness            0.353699051       0.9761171
## dorsolateralprefrontal.rh.area           0.494877863       0.9761171
## occipital.rh.area                        0.754968477       0.9761171
## occipital.rh.thickness                   0.501924742       0.9761171
## motorpremotor.rh.area                    0.804276395       0.9437441
## ventralfrontal.lh.thickness              0.970430749       0.9761171
## precuneus.rh.area                        0.414311663       0.9761171
## occipital.lh.thickness                   0.313896344       0.9761171
## medialtemporal.rh.thickness              0.947472757       0.9761171
## parsopercularis.lh.area                  0.622475737       0.9761171
## parsopercularis.rh.area                  0.692160839       0.9761171
## precuneus.lh.area                        0.672030293       0.9761171
## superiortemporal.lh.area                 0.863420323       0.9761171
##                                    fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## middletemporal.rh.thickness              0.4294997         0.8010153
## posterolateraltemporal.rh.area           0.2496220         0.2062780
## anteromedialtemporal.lh.area             0.2934903         0.2062780
## superiorparietal.rh.area                 0.4294997         0.4795132
## medialprefrontal.lh.thickness            0.3287141         0.8010153
## middletemporal.lh.thickness              0.2934903         0.8010153
## superiorparietal.lh.area                 0.4294997         0.4795132
## dorsomedialfrontal.lh.area               0.8575763         0.4795132
## medialprefrontal.rh.thickness            0.6237181         0.8010153
## temporalpole.rh.thickness                0.5091022         0.8704098
## orbitofrontal.rh.area                    0.1511874         0.8977969
## temporalpole.lh.thickness                0.3287141         0.4795132
## occipital.lh.area                        0.7988236         0.8010153
## posterolateraltemporal.lh.area           0.6237181         0.8010153
## inferiorparietal.lh.area                 0.2934903         0.8010153
## dorsolateralprefrotal.lh.thickness       0.5091022         0.8318132
## inferiorparietal.lh.thickness            0.6237181         0.4795132
## orbitofrontal.lh.area                    0.2261009         0.5694430
## perisylvian.lh.thickness                 0.6237181         0.9158838
## superiortemporal.rh.area                 0.8798973         0.8977969
## dorsolateralprefrotal.rh.thickness       0.9002974         0.8010153
## anteromedialtemporal.rh.area             0.6237181         0.8704098
## dorsomedialfrontal.rh.area               0.2496220         0.8010153
## inferiorparietal.rh.area                 0.9002974         0.8010153
## motor_premotor_SMA.rh.thickness          0.2496220         0.8010153
## motor_premotor_SMA.lh.thickness          0.1511874         0.8010153
## perisylvian.rh.thickness                 0.6237181         0.8010153
## ventralfrontal.rh.thickness              0.8798973         0.8704098
## ventromedialoccipital.lh.thickness       0.2261009         0.8760283
## inferiorparietal.rh.thickness            0.8907920         0.8010153
## ventromedialoccipital.rh.thickness       0.6237181         0.9158838
## superiorparietal.rh.thickness            0.9487756         0.8010153
## medialtemporal.lh.thickness              0.6237181         0.8010153
## dorsolateralprefrontal.lh.area           0.6237181         0.8010153
## motorpremotor.lh.area                    0.6237181         0.8010153
## superiorparietal.lh.thickness            0.7988236         0.8010153
## dorsolateralprefrontal.rh.area           0.8474733         0.8010153
## occipital.rh.area                        0.7988236         0.8977969
## occipital.rh.thickness                   0.7988236         0.8010153
## motorpremotor.rh.area                    0.6237181         0.8977969
## ventralfrontal.lh.thickness              0.8798973         0.9704307
## precuneus.rh.area                        0.7988236         0.8010153
## occipital.lh.thickness                   0.8798973         0.8010153
## medialtemporal.rh.thickness              0.8798973         0.9676318
## parsopercularis.lh.area                  0.6237181         0.8704098
## parsopercularis.rh.area                  0.8798973         0.8743084
## precuneus.lh.area                        0.7988236         0.8718231
## superiortemporal.lh.area                 0.6237181         0.9158838

Make effect size maps

dotSize = 6
fontSize=28

fdr_thresh = 0.05
# output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
# output_res_fdr[order(output_res_fdr$pval),]

#------------------------------------------------------------------------------
# df2plot = output_res[4:nrow(output_res),]
df2plot = output_res
df2plot$region = rownames(df2plot)
new_region_labels = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal")
hemi_labels = c("LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH")
metric_labels = c("SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT")

df2plot$region = new_region_labels
df2plot$hemisphere = hemi_labels
df2plot$metric = metric_labels

colors2use = get_ggColorHue(7)

#------------------------------------------------------------------------------
# Surface Area TD vs ASD Good Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.TD_vs_Good),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.TD_vs_Good),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "SA_TD_vs_ASDGood_effectSizeMaps.pdf"), plot = pf)
pf

#------------------------------------------------------------------------------
# Surface Area TD vs ASD Poor Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.TD_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.TD_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "SA_TD_vs_ASDPoor_effectSizeMaps.pdf"), plot = pf)
pf

#------------------------------------------------------------------------------
# Surface Area ASD Good vs ASD Poor Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.Good_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="SA")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenAr,
          mapping=aes(fill=es.Good_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "SA_ASDGood_vs_ASDPoor_effectSizeMaps.pdf"), plot = pf)
pf

#------------------------------------------------------------------------------
# Cortical Thickness TD vs ASD Good Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.TD_vs_Good),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.TD_vs_Good),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "CT_TD_vs_ASDGood_effectSizeMaps.pdf"), plot = pf)
pf

#------------------------------------------------------------------------------
# Cortical Thickness TD vs ASD Poor Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.TD_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.TD_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "CT_TD_vs_ASDPoor_effectSizeMaps.pdf"), plot = pf)
pf

#------------------------------------------------------------------------------
# Cortical Thickness ASD Good vs ASD Poor Effect size maps
# LH
# df2use = subset(d2plot, d2plot$hemisphere=="LH")
df2use = subset(df2plot, df2plot$hemisphere=="LH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pl = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.Good_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="left") +
  # scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

# RH
# df2use = subset(d2plot, d2plot$hemisphere=="RH")
df2use = subset(df2plot, df2plot$hemisphere=="RH" & df2plot$metric=="CT")
df2use$label2fill = factor(c(1:dim(df2use)[1]))

pr = ggseg(.data = df2use,
          atlas = chenTh,
          mapping=aes(fill=es.Good_vs_Poor),
          # mapping=aes(fill=Fstat),
          # mapping=aes(fill=label2fill),
          position="stacked",
          colour="black",
          hemisphere="right") +
  # scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) +
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits = c(-0.7,0.7)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  guides(fill=guide_legend(title="Effect Size")) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
ggsave(filename = file.path(plotdir, "CT_ASDGood_vs_ASDPoor_effectSizeMaps.pdf"), plot = pf)
pf

Make A-P and D-V and Effect Size Gradient Plots

regions2use = colnames(fs_data)[5:(ncol(fs_data)-4)]
full_regnames = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal")

sa_ordering = c(1,12,7,10,
                4,6,9,3,
                8,11,2,5)
ct_ordering = c(1,7,6,3,
                10,4,11,5,
                8,2,12,9)
gensim_order = c(sa_ordering,sa_ordering, ct_ordering,ct_ordering)

df4plot = output_res[regions2use,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")]
df4plot$phenotype = "CT"
df4plot$phenotype[1:24] = "SA"
df4plot$phenotype = factor(df4plot$phenotype) 

df4plot$hemisphere = "RH"
df4plot$hemisphere[1:12] = "LH"
df4plot$hemisphere[25:36] = "LH"
df4plot$hemisphere = factor(df4plot$hemisphere) 

df4plot$full_region_names = factor(full_regnames)
df4plot$gensim_gradient = gensim_order

ct_df = subset(df4plot, df4plot$phenotype=="CT")
sa_df = subset(df4plot, df4plot$phenotype=="SA")

ct_df$cluster = "dorsal"
ct_df$cluster[c(6,7,8,9,10,11,12)] = "ventral"
ct_df$cluster = factor(ct_df$cluster)

sa_df$cluster = "anterior"
sa_df$cluster[c(6,7,8,9,10,11,12)] = "posterior"
sa_df$cluster = factor(sa_df$cluster)

dotSize = 10
fontSize = 18

p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Good
## t = -1.5876, df = 22, p-value = 0.1267
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6410888  0.0950901
## sample estimates:
##        cor 
## -0.3206047
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 0.026067, df = 21.915, p-value = 0.9794
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1712679  0.1756271
## sample estimates:
##  mean in group dorsal mean in group ventral 
##          -0.007408816          -0.009588433
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Poor
## t = -0.79935, df = 22, p-value = 0.4326
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5351309  0.2525090
## sample estimates:
##        cor 
## -0.1680004
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = -1.0499, df = 12.781, p-value = 0.3132
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4389992  0.1521861
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           -0.05882051            0.08458602
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.Good_vs_Poor
## t = 0.66563, df = 22, p-value = 0.5126
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2786875  0.5147273
## sample estimates:
##      cor 
## 0.140505
p = ggplot(data = ct_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = -1.2707, df = 10.307, p-value = 0.2318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.36781016  0.09995858
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           -0.04447174            0.08945405
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Good
## t = 0.054893, df = 22, p-value = 0.9567
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3935522  0.4131487
## sample estimates:
##        cor 
## 0.01170234
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 1.1522, df = 16.828, p-value = 0.2654
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08611438  0.29297353
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.03023708             -0.07319250
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Poor
## t = -0.57123, df = 22, p-value = 0.5736
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4999111  0.2969860
## sample estimates:
##        cor 
## -0.1208943
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = 1.6048, df = 11.985, p-value = 0.1345
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07913237  0.52140838
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.09699741             -0.12414060
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.Good_vs_Poor
## t = -0.58651, df = 22, p-value = 0.5635
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5023311  0.2940367
## sample estimates:
##        cor 
## -0.1240773
p = ggplot(data = sa_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = 0.71299, df = 9.2556, p-value = 0.4934
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2262792  0.4358494
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.06145044             -0.04333464